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SUMMARY 

The  governing  equations  for  compressible  three-dimensional  turbulent 
boundary  layers  in  a  general  coordinate  system  are  given.  Methods  for  the 
solution  of  these  equations  and  their  integral  counterparts  are  described 
and  compared.  The  available  finite  difference  techniques  are  discussed  in 
detail . 
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INTRODUCTION 


1  would  Like  to  begin  by  saying  what  an  honour  and  a  pleasure  it  is  i or  me  u>  ; i v> 

this  review  lecture.  To  some  extent  my  pleasure  is  enhanced  by  a  sense  of  d£ja  vu  as  it 

is  almost  exactly  ten  years  since  I  spoke  on  a  similar  topic  at  this  very  Ir.stituti. 

The  occasion  was  then  EUROMECH  33  '  on  three-dimensional  turbulent  boundary  layers.  At 

that  meeting  I  described  the  integral  calculation  method^  I  had  developed  lor  compress:!,  i e 

turbulent  boundary  layers  and  laid  particular  stress  upon  the  generality  of  Lhe  coordinate 

system  I  had  adopted.  At  that  time  the  method  was  probably  the  most  advanced  (at  least  i r. 

its  range  of  applicability)  available  and  will  serve  as  a  useful  datum  against  which  t. 

measure  progress  during  the  past  decade.  Of  course,  as  I  hope  to  show,  real  progress  has 

been  made.  Other  integral  methods  have  adopted  the  general  coordinate  system  and  have 

used  more  sophisticated  assumptions  about  the  velocity  profiles  and  entrainment.  The 

major  advance  has  been,  however,  in  the  development  of  the  finite-difference  methods  and 

I  hope  to  spend  some  time  describing  these  in  detail.  As  a  guard  against  complacency  let 

me  state  here  however,  that  progress  has  not  been  achieved  in  the  types  of  turbulence 

3 

modelling  we  use  in  these  methods.  In  1972  Wesseling  and  Lindhout  descriV.  .  their 

finite-difference  method  which  although  confined  to  incompressible  flow  an a  Cartesian 
.  4 

coordinates  used  Bradshaw's  vector  equation  for  the  turbulent  shear  stress.  Ten  years 
later  almost  all  the  finite-difference  methods  appear  to  have  reverted  to  the 
primitive  eddy-viscosity  concept  of  turbulence  modelling. 


EUROMECH  33  also  marked  the  birth  of  another  significant  development  in  this  field, 
the  organisation  of  a  series  of  meetings  to  compare  results  from  the  various  calculations 
methods  with  each  other  and  where  possible  against  experiment.  The  first  of  these 
meetings,  EUROMECH  60,  became  known  as  the  Trondheim  Trials"’  and  prompted  the  EUROVISC 
Committee  of  national  representatives  to  form  a  working  party  which  has  since  organised 
two  further  meetings,  at  Stockholm**  in  1978  and  Amsterdam''  in  1979.  Finally,  as  you  will 
all  know  a  third  meeting  will  immediately  follow  this  Symposium  so  that  any  predictions 
I  may  be  foolish  enough  to  make  concerning  the  capabilities  ol  the  various  calculation 
methods  will  rapidly  be  put  to  the  test  and  almost  certainly  shown  to  be  false! 


As  with  any  review  paper  little  of  what  follows  is  original.  In  particular  I  have 
drawn  heavily  upon  the  works  of  J.P.F.  Lindhout''  of  the  SLR  and  my  colleague  at  the  RAE, 
Caroline  Boyd  .  To  them  and  also  to  all  those  whose  names  appear  in  the  list  of 
references  I  offer  my  sincere  thanks. 

2  THE  THREE-DIMENSIONAL  BOUNDARY  LAYER  EQUATIONS0*" '  1 

For  general  appl icab i  1 i ty  most  three-dimensional  boundary— layer  calculation  ruthods 
use  a  non-orthogonal  curvilinear  coordinate  system;  x,  y  in  the  surface,  r.  normal  to 
the  surface  with  an  angle  A  (x,y)  between  the  x  and  y  directions.  In  this 
coordinate  system  the  continuity,  momentum  and  energy  equations  may  be  written  as. 
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Here 


pw  =  pw  +  p  V 


An  element  of  length  on  the  surface  ds  is  given  by 

ds^  =  h^dx^  +  h^dy^  +  2h)h2  cos  X  dxdy  (5) 

and  the  total  velocity  U  is  given  by 

=  u^  +  v^  +  2uv  cos  X  (6) 

where  the  metric  coefficients  h|  and  h9  are  to  the  boundary  layer  approximation 
functions  of  x  and  y  but  not  z  ,  ie 


h, 

=  h j (x,y)  , 

h2  =  h2(x,y)  . 

The  boundary  conditions  are 

z  =  0  ;  u  - 

v  =  0  ,  w 
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The  curvature  parameters  kj,  k^,  k^  and  k^  are  8iven  by 
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For  an  orthogonal  system,  A  =  90°,  kj  =  -k)2,  k2  =  -k2 1  . 

Rather  than  solve  the  energy  equation  (4),  for  adiabatic  flow  the  density  at  each 
point  in  the  boundary  layer  is  frequently  assumed  to  be  approximated  by  the  Crocco 
relation 


_ \±(h  ,,_V| 

h^  sin  X  |^3x  2  C0S  A  3y  J 

[w  (hi  cos  x)  "  "5t] 


hjh7  sin  X 


‘12  s 


1 

C21  sin  X 


itt  [-  (ki  *  57  f;)  *  005  (k>  *  r;  ij)] 

[-  (k2  c“ 


pe 

c 


( !  2 ) 


where  r  ,  the  recovery  factor  is  taken  as  0.84  or  0.85  for  laminar  flows  and  0.89  for 
turbulent  flow. 

At  the  boundary-layer  edge  equations  (2)  and  (3)  reduce  to 
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The  two  Euler  equations,  (13)  and  (14),  may  be  used  to  eliminate  p  from  equations  i  ) 
and  (3),  or  for  a  given  pressure  distribution  p  (x,y)  may  be  solved  for  ug  and  v 
subject  to  appropriate  initial  and  boundary  conditions. 

3  INTEGRAL  METHODS 

The  earliest  and  simplest  attempts  to  find  solutions  to  the  general  equations  (1)  to 

(3)  involved  the  use  of  integral  methods  in  which  the  equations  are  integrated  term  by 

term  across  the  boundary  layer  from  the  wall  to  the  outer  edge.  These  integral  equations 

•  >2 

were  first  derived  by  Myring  .  To  proceed  further  it  is  then  necessary  to  assume  some 
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form  for  the  velocity  profiles  in  the  boundary  layer.  All  integral  methods  make  use  of 

streamwise  and  crosswise  profiles.  That  is  profiles  in  the  directions  of  the  external 

streamlines  (s)  and  their  orthogonal  trajectories  (n)  respectively.  It  was  therefore 

necessary  for  Myring  to  derive  the  transformation  relationships  between  the  integral 

thicknesses  in  the  x,  y  coordinate  system  and  the  corresponding  thicknesses  in  the  s,  n 

coordinate  system.  The  equations  may  however  be  derived  rather  more  directly  by  starting 
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from  the  integral  equations  in  streamline  coordinates  and  transforming  the  independent 
variables  s,  n  to  x,  y  as  shown  below.  This  derivation  will  be  useful  when  we  come 
to  discuss  inverse  methods. 

The  integral  equations  in  streamline  coordinates  are: 
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where  C,,  is  the  entrainment  coefficient. 

E 

These  equations  are  transformed  to  the  general  x,  y  non-or thogona 1  coordinate  svstem  bv 
using  the  operators 
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where  >  is  the  angle  between  the  x  and  y  directions  and  a  the  angle  between  the 
x  and  s  directions  so  that 
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The  curvature  terms  k  and  k  are  derived  as  follows.  For  k  we  equate  the 

s  n  s 

continuity  equation  (I)  for  the  external  flow  in  the  x,  y  coordinate  system  to  its 
counterpart  in  the  s,  n  coordinate  system.  The  result  is 


k  =  -  — r — -  [—  f  h  s  i  n  ( '•  -  x )]  +  - —  (h  i.in  x)l  .  <  -  1  ) 

s  l^h,,  sin  *  L  'x  1  2  j  ,y  1  J 
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with  the  result:  thaL 
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We  note  from  equation  (.23)  that  if  the  extern.il  flow  is  irrot.u ional  then 

k  =  tl/L'  )  OF  /  ~>n)  and  some  simplification  of  equations  (15)  to  (17)  results, 
nee  n 

Up  to  this  point  no  approximations  have  been  made  but  to  proceed  further  it  is 
necessary  to  assume  forms  for  the  velocity  profiles  in  order  to  derive  relationships 
between  t1 e  unknowns  so  that  they  may  be  reduced  in  number  to  three.  The  various  assump¬ 
tions  which  have  been  made  are  given  below. 

2  .  .  l  s 

Smith  assumed  power  law  pro: lies  for  the  streamwiso  flow  and  Mager  '  profiles  for 

the  crossflow.  Streamwiso  skin  friction  and  the  entrainment  coefficient  were  derived 

from  auxiliary  relationships.  Later  the  method  was  modified  so  that  the  entrainment 

coefficient  was  calculated  bv  an  extension  of  the  lag-entrainment  method  of 
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Green,  <■ 

Stock  ~  modified  Smith's  method  so  that  the  streamwiso  profile  was  represented  by 
Coles  law  of  the  wall  plus  wake.  This  removed  the  need  for  an  empirical  skin  friction 
relationship.  In  addition  the  entrainment  was  calculated  from  Horton's  lag  entrainment 
relationship. 

18 

Cross  uses  a  three-dimensional  form  of  the  law  of  the  wall  and  wake  together  with 
Thompson*  s  entrainment  relationship. 


r 
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Cous teix * ^ and  Aupoix  use  velocity  profiles  derived  from  an  analysis  of  similarity 
solutions . 

The  methods  are  listed  above  in  order  of  increasing  complexity.  They  all  use  very 
similar  numerical  methods  to  solve  the  equations  and  it  is  interesting  to  note  that  this 
increase  in  complexity  is  reflected  in  increased  computing  times  as  listed  by  Lindhout 
for  the  Amsterdam  Workshop  test  case.  The  times  are  shown  in  Table  I  in  terms  rf  effort 
measured  in  Mega  Floating  Point  Operations  (.MFLOP)  per  surface  point. 

Table  I 


Method 

r 

Effort  per  point  (MFLOP) 

Smith 

0.007 

Stock 

0.013 

C  ross 

0.018 

Cousteix  and  Aupoix 

0.033 

These  figures  must  be  interpreted  with  caution.  The  methods  have  usually  been  written 
with  clarity  and  ease  of  alteration  rather  than  speed  as  the  prime  objective  and  in  addi¬ 
tion  the  numerical  method  of  Cousteix  and  Aupoix  is  more  complicated  chan  those  by  the 
others.  The  important  point  to  note  is  that  all  the  methods  are  so  fast  that  realistic 
configurations  can  be  calculated  in  very  acceptable  run  times  and  even  the  slowest 
integral  method  remains  an  order  of  magnitude  faster  than  a  finite  difference  method. 

The  assumption  of  velocity  profiles  and  where  necessary  additional  empirical 
relationships  reduces  the  numbei  of  unknowns  to  three  (A,  B,  C  say)  and,  typically,  three 
equations  of  the  form 
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where  all  the  boundary  layer  parameters  are  now  known  functions  of  A,  B,  C  and  U  . 

1  2 

The  equations  are  hyperbolic  with  three  real  characteristics  and  it  is  necessary 
to  take  account  of  these  characteristic  directions  when  evaluating  the  y  derivatives  as 
shown  in  Fig  1.  To  a  good  approximation  the  upper  and  lower  bounds  of  the  characteristic 
directions  with  respect  to  the  x  axis  are  given  by  a  and  a  +  (5  respectively  where  e 
is  the  direction  of  the  wall  streamline  with  respect  to  s  .  With  the  y  derivatives 
evaluated  in  this  manner  the  solution  is  advanced  in  the  x  direction  by  a  simple 
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explicit  method.  Smith"  for  example  uses  a  two  step  Euler  method  whilst  Cousteix  and 

Aupoix'^’^  use  a  fourth  order  Runge  Kutta  method.  The  explicit  nature  of  the  method 

restricts  the  forward  step  size  to  that  allowed  by  the  Courant  Friedrichs  Levy  condition 
but  as  noted  above  run  times  are  short  and  no  attempt  has  been  made  to  remove  this  res¬ 
triction  by  developing  an  implicit  method. 

4  DISPLACEMENT  THICKNESS  AND  EQUIVALENT  SOURCE  DISTRIBUTION 

For  calculations  involving  viscous-invisci d  interactions  the  effect  of  the  boundary 
layer  upon  the  external  flow  may  be  represented  as  a  change  in  boundary  condition  for  the 
inviscid  calculation  from  one  of  zero  velocity  normal  to  the  surface  to  one  in  which  the 

velocity  component  normal  to  the  surface,  m  ,  is  given  by 
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m  is  often  called  the  transpiration  velocity  or  equivalent  source  strength. 


Alternatively  the  effect  of  the  boundary  layer  may  be  represented  as  a  displacement 


of  the  surface  through  a  distance  6*  normal  to  the  surface  where 


is  given  by 
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an  additional  partial  differential  equation  which  must  be  solved  if  this  approach  is 
adopted . 

5  INVERSE  INTEGRAL  METHODS 

For  the  calculation  of  separated  flows  it  is  necessary  to  invert  the  method  so  that 
the  boundary  layer  development  is  specified  and  the  velocity  distribution  is  calculated. 
We  can  readily  see  that  if  two  of  the  boundary  layer  thicknesses  (B,  C  say)  are  specified 
as  functions  of  x  and  y  in  equations  (25)  then  the  equations  may  be  rearranged  to  the 
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Solutions  of  this  type  for  infinite  yawed  wings,  in  which  (3/3y)  =  0  ,  have  been 
2122 

published  by  Cousteix  and  Stock  .  In  general,  however,  the  inverse  method  will  be  part 
of  an  interactive  calculation  and  only  the  source  strength  m  or  the  displacement  thick¬ 
ness  6*  will  be  available  as  input  to  the  inverse  method.  An  additional  equation  is 


then  required  and  this  is  given  by  equation  (22).  The  inverse  method  is  then  a  simul¬ 
taneous  solution  of  equations  (25),  (26)  or  (27)  and  (22).  Whilst  this  is  obviously 
possible  in  principle  the  characteristic  directions  may  cause  difficulties  for  any  march¬ 
ing  method  and  no  working  method  of  this  type  has  yet  been  published. 

6  FINITE-DIFFERENCE  METHODS 


The  alternative  to  the  integral  approach  is  to  approximate  the  terms  of  equations 

(I)  to  (4)  by  finite  differences.  If  this  is  done  in  a  consistent  manner  the  solution  of 

the  difference  equations  will  hopefully  converge  to  the  solution  of  the  differential 

equations  as  the  step  sizes  are  reduced  to  zero.  The  earliest  methods  strongly  resembled 

the  numerical  techniques  adopted  for  integral  methods  in  that  they  were  explicit.  That  is 

to  obtain  the  solution  at  the  point  x.,  v.,  z,  the  solution  is  written  in  terms  of  the 

l  '  j  k 

unknowns  at  x.,  v.,  z,  and  known  values  in,  the  y,  z  plane  at  x.  ,  .  Fpwind 

l  ‘  J  k  l-l 

differencing  is  usually  used  in  the  v  direction  as  shown  in  Fig  1.  More  recently  the 
trend  has  been  to  implicit  methods  in  which  the  solution  is  sought  simultaneously  along  a 
complete  z  column  at  x^,  y ^  .  Once  again  upwind  differencing  is  used  in  the  y  direc¬ 
tion  and  central  differences  are  used  to  evaluate  the  second  derivative  which  occurs  in 

25 

the  z  direction.  Typical  of  this  type  of  method  is  that  of  Chang  and  Patel  .  It  will 

be  seen  from  Fig  1  that  no  stability  restrictions  arise  from  the  CFI.  condition  for  this 

approach  but  the  method  is  restricted  to  flows  in  which  v  ,  the  velocity  in  the  v 

26 

direction,  is  everywhere  positive.  If  v  is  positive  the  methods  of  l.indhout  and  Nash 
27 

and  Scruggs  use  differencing  techniques  identical  to  those  of  Chang  and  Patel.  Those 
28 

used  by  Mclean  are  similar  but  in  addition  Mclean  uses  three  point  backward  differences 
in  an  attempt  to  preserve  second  order  accuracy.  Once  v  becomes  negative  it  is  necessary 

to  evaluate  the  y  derivative  at  the  k.  .  station  and  then  these  methods  are  subject  to 

l-l  J 

the  stability  restriction  that  the  characteristic  (streamline)  through  x.,  y.,  z  must 

■  n  t  1  k 


intersect  x^_,  ,  j 


jqt  'J  k 

y  -  j  +1  .  The  zig-zag  technique  used  by  Kordulla"',  and  originally 

30  •  _  , _  _ _ _  •  _  ..  . _ 


formulated  by  Krause  ,  is  similar  but  here  rather  than  switch  operators  as  v  changes 
sign  the  same  second  order  operator  is  used  throughout  as  shown  in  Fig  2.  In  addition 
second  order  accuracy  in  the  x  direction  is  achieved  by  centering  the  difference 


equations  on  M  , 

(xi-i’yj’V 

may  be  written 

/  3  (  ) 
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If  varying  step  sizes  Ay  j ,  hy^  say  are  used  then  second  order  accuracy  is  only  preserved 
if  Ay^-Ayj  is  kept  of  order  Ay ^  •  Ay9  .  Similar  restrictions  apply  to  varying  step 

sizes  in  the  z  direction.  In  order  to  overcome  these  restrictions  and  permit  arbitrarily 

.  •  •  .  3 1 

varying  step  sizes  whilst  preserving  second  order  accuracy  Cebeci  has  used  the  box 

scheme.  Here  as  shown  in  Fig  3  the  equations  are  centred  on  M(x.  ,,v.  ,,z  ,)  and 

1  i  '  J _  1  k  * 

(  ).  t  .  j  .  I  is  taken  as  the  average  of  the  eight  corner  values.  Derivatives  are 
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written  as 
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The  second  derivative  in  the  z  direction  is  obtained  by  defining  new  unknowns  f,  g 

i  7  2 

say  equal  to  3u/dz  and  3v/dz  .  Then  3“ u/3z“  =  7«€/3z  and  3~v/2 z~  =  -g  /  'z  . 

The  box  technique  can  only  be  used  if  v  is  positive  and  originally  Cebeci,  •  J 

proposed  a  solution  technique  in  which  the  solution  marched  in  either  the  positive  or 

negative  y  directions  according  to  the  sign  of  v  .  This  technique  poses  programming 
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problems  (which  appear  to  have  been  overcome  in  the  method  of  Lindhout  to  which  wo  will 

,  P 

shortly  return)  and  so  Cebeci,  •  :  a.  ~  proposed  two  alternative  techniques.  The  first  of 

these,  the  zig-zag  box  scheme,  combines  the  zig-zag  and  box  methods  when  v  is  negative 

so  that  M  becomes  x.  . ,  y. ,  z,  ,  and  now  derivatives  are  written  as 

1—4  i  k-J 


fifil  =  ^1-  f(  ).  . 
L  ^  JM  2“xi-i  l 


J,k  +  (  )i,j,k-l  (  }i-l ,j,k  (  } 


i-l  .  j  .k-l] 


•  t[{<  h,j.k*  <  >i.j,k-i  -  *  > 

+  {(  }  i-l ,j+l ,k  +  (  ) i-l ,j  +  l ,k 


k  +  (  } i  ,  j ,k- 1  "  (  )i,j-l,k  (  } i , j- 1 ,k- 1 }/"yj-l 


k-I  ^  ^  ^ i- 1 , j ,k- 1 f / j i 


and 
(  ). 


[  3z  jM  "  2AzR_|  [l  }i,j,k  +  (  )i-l,j,k  "  (  } i , j , k- 1  ‘  )i-l,j,k-l] 

4  ^ i , j ,k  *•  ^i,j,k-l  +  ^  ^i-l,j,k  +  1  'i-l,j,k-l] 


The  use  of  these  approximations  does  however  re-introduce  the  stability  restriction 
mentioned  above  together  with  the  restriction  on  the  variation  of  the  y  step  sizes  if 
second  order  accuracy  is  to  be  preserved.  Cebeci,  r*  ,7.'  other  method  in  Ref  32  is  the 


ill  the  momentum 


characteristic  difference  scheme.  Fig  4,  in  which  the  inertia  terms 
equations  are  written  as 


4x 


where  2  is  the  direction  of  the  streamline  or  characteristic  through  x,  y,  . 

If  P.  .  ,  |  cuts  x._.  constant,  z  .  constant,  at  Q  the  equations  centred  now  on 

PQ/2  where  P  =  x.,  x.,  z,  .  can  be  differenced  as  ior  example 
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This  scheme  cannot  however  be  simply  applied  to  the  continuity  equation  which  is  usually 
dealt  with  by  a  zig-zag  scheme.  For  the  characteristic  scheme  to  remain  second  order 
accurate  i  must  be  evaluated  at  PQ/2  and  u  must  be  accurate  to  within  0  f 3! . 

q 

This  last  condition  can  be  satisfied  by  using  quadratic  interpolation  through  three 

Q 

points  say  y j _ ]  ,  y^  and  Yj  +  i  •  Boyd  has  compared  these  three  C.ebeci  difference 
schemes  for  the  Trondheim  Trials^  and  NLR  test  cases  and  finds  the  characteristic  scheme 
to  be  less  robust  and  to  require  roughly  three  times  the  computation  time.  Boyd's 
preferred  strategy  is  to  use  the  box  scheme  unless  any  point  on  a  z  column  requires  th< 
zig-zag  in  which  case  the  zig-zag  scheme  is  used  for  all  that  column.  Cebeci,  • 
found  that  the  characteristic  scheme  allowed  them  to  find  solutions  closer  to  separation 
than  the  other  schemes  but  Boyd  found  that  for  the  NLR  test  case  these  differences  were 
very  slight. 

It  was  mentioned  above  that  Cebeci  originally  adopted  a  marching  scheme  which 
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depended  upon  the  sign  of  v  .  Lindhout  uses  just  such  a  scheme  and  has  devised  an 
algorithm  to  define  the  optimal  marching  scheme.  The  advantage  of  this  scheme  is  that 
the  stability  restriction  on  the  x  step  is  minimised,  since  the  explicit  operators 
tend  only  to  be  used  as  v  changes  sign  and  is  therefore  close  to  zero. 


The  end  result  of  any  of  the  implicit  methods  is  a  set  of  non-linear  algebraii 
equations  which  are  linearised  by  the  use  of  iteration  and  then  solved  by  standard  rnrtr:-. 
techniques.  The  need  to  use  iteration  makes  comparisons  of  timings  between  the  metnous 

particularly  difficult.  Some  comparative  timings  are  given  in  Ref  7  but  they  are  not 

8 

reproduced  here  as  Boyd  has  found  that  the  number  of  x  steps  in  one  example  could  be 
doubled  with  Htde  or  no  increase  in  run  times.  This  apparent  anomaly  occurs  because 
the  smaller  step  size  produces  a  better  initial  guess  at  the  next  x  station  and  hence 
fewer  iterations  are  required.  A  similar  mechanism  may  explain  why  Lenmerman  found 
that  even  the  old  explicit  methods  remained  competitive. 

The  stability  restriction  may  be  removed  altogether  if  the  method  is  made  implicit 

in  both  the  y  and  z  directions  as  is  done  by  Nash  and  Scruggs"^,  Patankar  and 

34  35  ... 

Spalding  and  Baker  .  A  simultaneous  solution  technique  is  then  required  for  the  whole 

y,  z  plane.  Nash  and  Scruggs  use  an  alternative  direction  implicit  method,  :  aakar  and 

Spalding  their  simple  algorithm  and  Baker  a  finite-element  approach.  The  use  oi  such  a 

fully  implicit  technique  allows  the  inclusion  of  the  diffusion  in  the  y  direction  and 

this  has  been  done  by  both  Baker  and  Spalding.  In  addition  if  one  includes  the  z-di recti 

momentum  equation  the  pressure  may  be  allowed  to  vary  with  z  although  i;  becomes 

necessary  to  uncoupLe  the  pressure  in  the  x  momentum  equation  from  that  in  the  y  and 

z  momentum  equations.  For  further  details  see  Ref  34. 

Comparisons  between  the  results  of  the  various  difference  methods  may  he  found  in 
Refs  5,  6,  7,  8  and  33  but  as  mentioned  above  all  timings  contained  therein  must  be  treat 
with  extreme  caution  and  certainly  no  clear  'best  buy'  has  emerged  to  date. 

g 

No  inverse  finite-difference  method  has  been  published  but  Boyd  has  suggested  how 
a  method  using  equations  (1)  to  (4)  or  (12),  (22)  and  (26)  could  be  devised. 

7  CONCLUSIONS 

Calculation  methods  for  three-dimensional  turbulent  boundary  layers  on  practical 
shapes  are  now  fairly  widely  available.  Of  the  two  types,  integral  and  finite-difference 
the  integral  methods  retain  their  advantages  of  speed  and  robustness  which  make  then 
particularly  suitable  for  viscous-inviscid  interaction  calculations.  The  finite- 
difference  methods  are  however  not  restricted  by  any  velocity  profile  assumptions  and 
allow  more  advanced  turbulence  models  to  be  more  readily  incorporated.  In  addition  thev 
hold  the  promise  that  they  rather  than  experiment  may  suggest  new  profile  families  for 
integral  methods.  The  emergence  of  inverse  methods  both  integral  and  finite-difference 
would  appear  to  be  imminent. 

Note  added  in  proof 

No  sooner  were  the  above  words  committed  to  print  than  an  inverse  finite- 

36 

difference  method  appeared.  This  method  ,  by  Formery  and  Delery,  like  the  integral 
methods  described  in  section  5,  requires  two  boundary  layer  thicknesses  to  be  specified. 
This  may  make  its  application  to  viscous  inviscid  interactive  calculations  difficult.  It 
is  not  however  restricted  to  infinite  yawed  wing  flows  and  good  agreement  is  shown  with 
both  infinite  yawed  wing  and  fully  three-dimensional  experiments.  Difficulties  produced 
by  the  characteristics  coming  from  downstream  are  avoided  in  this  method  by  using  the 
FLARE  approximation  in  which  for  example  the  term  ^u /(tx  is  set  to  zero  if  u  is  negati 
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